Engineering a computational surface modeling toolbox in Julia

Ray M. Yang

Washington University in St. Louis

2024-11-04

@ PChem Group Meeting

Overview

  1. Introduction to FreeBird.jl
    Design, functionality, and extensibility

  2. Introduction to Research Software Engineering (RSE)
    What is RSE, and how does it fit in Chemistry

What is FreeBird.jl

  • It’s a computational package
  • It’s developed by the Wexler Group
  • It does surface free energy calculations via nested sampling
  • It’s written in Julia, a modern, high-level programming language

What is Julia?

is a high-level, high-performance, dynamic programming language for technical computing. It’s fast, easy to learn, and has a syntax that is familiar to users of other technical computing environments.

  • Version 1.0 was released in 2018
  • Developed by mathematician Alan Edelman et al. at MIT
  • It’s designed for high-performance numerical analysis and computational science
  • It’s free and open source

FreeBird.jl design

  • It’s modular
  • It’s extensible
  • It’s user-friendly
  • It’s well-documented
  • It’s open-source

FreeBird.jl functionality

  • Systems: atomistic and lattice models
  • Energy calculator: Lennard-Jones and lattice Hamiltonians
  • Methods: nested sampling, Metropolis Monte Carlo, exact enumeration, Wang-Landau sampling

Nested sampling

Atomistic walkers

using FreeBird
walkers = AtomWalker.(generate_initial_configs(120, 562.5, 6))
lj = LJParameters(epsilon=0.1, sigma=2.5)
LJParameters(0.1 eV, 2.5 Å, Inf, 0.0 eV)

Atomistic walkers

walkers[1]
AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

Atomistic walkers

walkers[1].configuration
FastSystem(H₆, periodic = FFF):
    bounding_box      : [      15        0        0;
                                0       15        0;
                                0        0       15]u"Å"

    AtomView(H,  [ 4.01502,  4.18315,  10.6918]u"Å")
    AtomView(H,  [ 4.94146,  3.03619,  4.64721]u"Å")
    AtomView(H,  [ 5.35381,  11.0872,  7.85456]u"Å")
    AtomView(H,  [ 9.05429,  1.99427,  7.42124]u"Å")
    AtomView(H,  [  13.382,  13.0819, 0.368069]u"Å")
    AtomView(H,  [ 13.3007,  11.1216, 0.182911]u"Å")

             .------------------------------------.  
            /|                                    |  
           / |                                    |  
          /  |                                    |  
         /   |                                    |  
        /    |                                    |  
       /     |                                    |  
      /      |                                    |  
     /       |                                    |  
    *        |                                    |  
    |        |                                    |  
    |        |          H                         |  
    |        |  H                                 |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |              H                     |  
    |        .------------------------------------.  
    |       /                                H   /   
    |      /      H                        H    /    
    |     /                                    /     
    |    /                                    /      
    |   /                                    /       
    |  /                                    /        
    | /                                    /         
    |/                                    /          
    *------------------------------------*           

Construct the liveset

ls = LJAtomWalkers(walkers, lj)
LJAtomWalkers(AtomWalker{1}, LJParameters):
[1] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : 5.2747930136178764 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[2] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.03877740224016672 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[3] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.04244190870170766 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[4] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.10913780161856197 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[5] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.11327569913084067 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

⋮
Omitted 110 walkers
⋮

[116] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.07630973295538208 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[117] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : 1.3074885636307922 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[118] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : 132.93929690640402 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[119] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : 31.242106865621327 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

[120] AtomWalker{1}(
    configuration      : FastSystem(H₆, periodic = FFF, bounding_box = [[14.999999999999998, 0.0, 0.0], [0.0, 14.999999999999998, 0.0], [0.0, 0.0, 14.999999999999998]]u"Å")
    energy             : -0.009495483209446563 eV
    iter               : 0
    list_num_par       : [6]
    frozen             : Bool[0]
    energy_frozen_part : 0.0 eV)

LJParameters(0.1 eV, 2.5 Å, Inf, 0.0 eV)

Run nested sampling

mc = MCRandomWalkClone()
ns_params = NestedSamplingParameters(200, 0.1, 0.01, 1e-5, 1.0, 0, 200)
save = SaveEveryN(n_traj=10, n_snap=10_000)
energies, liveset, _ = nested_sampling_loop!(ls, ns_params, 10_000, mc, save)

Nested sampling results

liveset.walkers[1].configuration
FastSystem(H₆, periodic = FFF):
    bounding_box      : [      15        0        0;
                                0       15        0;
                                0        0       15]u"Å"

    AtomView(H,  [ 2.92958,  6.58633,  3.60833]u"Å")
    AtomView(H,  [  2.1724,  6.20473,  0.94135]u"Å")
    AtomView(H,  [0.227731,  6.37336,  2.94226]u"Å")
    AtomView(H,  [ 3.21941,  8.66147,  1.76325]u"Å")
    AtomView(H,  [0.510936,  8.44946,  1.10138]u"Å")
    AtomView(H,  [   1.269,  8.83089,   3.7619]u"Å")

             .------------------------------------.  
            /|                                    |  
           / |                                    |  
          /  |                                    |  
         /   |                                    |  
        /    |                                    |  
       /     |                                    |  
      /      |                                    |  
     /       |                                    |  
    *        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |        |                                    |  
    |       H|                                    |  
    |        .------------------------------------.  
    |       /  H                                 /   
    |   H  /     H                              /    
    |     H                                    /     
    |    /   H                                /      
    |   /                                    /       
    |  /                                    /        
    | /                                    /         
    |/                                    /          
    *------------------------------------*           

Nested sampling results

Nested sampling results

Run nested sampling

Nested sampling in 7 lines of code!

ls = LJAtomWalkers(AtomWalker.(generate_initial_configs(120, 562.5, 6)), LJParameters())
mc = MCRandomWalkClone()
ns_params = NestedSamplingParameters(200, 0.1, 0.01, 1e-5, 1.0, 0, 200)
save = SaveEveryN(n_traj=10, n_snap=10_000)
energies, liveset, _ = nested_sampling_loop!(ls, ns_params, 10_000, mc, save)

Or, in a single line if you really want to!

energies, liveset, _ = nested_sampling_loop!(LJAtomWalkers(AtomWalker.(generate_initial_configs(100, 100.0, 6)), LJParameters()), NestedSamplingParameters(200, 0.1, 0.01, 1e-5, 1.0, 0, 200), 10_000, MCRandomWalkClone(), SaveEveryN(n_traj=10, n_snap=10_000))

Accessing the documentation

help?> nested_sampling_loop!
search: nested_sampling_loop! nested_sampling_step!

  nested_sampling_loop!(liveset::AtomWalkers, ns_params::NestedSamplingParameters, n_steps::Int64, mc_routine::MCRoutine; args...)

  Perform a nested sampling loop for a given number of steps.

  Arguments
  ≡≡≡≡≡≡≡≡≡

    •  liveset::AtomWalkers: The initial set of walkers.

    •  ns_params::NestedSamplingParameters: The parameters for nested sampling.

    •  n_steps::Int64: The number of steps to perform.

    •  mc_routine::MCRoutine: The Monte Carlo routine to use.

  Returns
  ≡≡≡≡≡≡≡

    •  df: A DataFrame containing the iteration number and maximum energy for each step.

    •  liveset: The updated set of walkers.

    •  ns_params: The updated nested sampling parameters.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  nested_sampling_loop!(liveset::LatticeGasWalkers, ns_params::LatticeNestedSamplingParameters, n_steps::Int64, mc_routine::MCRoutine; args...)

  Perform a nested sampling loop on a lattice gas system for a given number of steps.

  Arguments
  ≡≡≡≡≡≡≡≡≡

    •  liveset::LatticeGasWalkers: The initial set of walkers.

    •  ns_params::LatticeNestedSamplingParameters: The parameters for nested sampling.

    •  n_steps::Int64: The number of steps to perform.

    •  mc_routine::MCRoutine: The Monte Carlo routine to use.

  Keyword Arguments
  ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡

    •  args...: Additional arguments.

  Returns
  ≡≡≡≡≡≡≡

    •  df: A DataFrame containing the iteration number and maximum energy for each step.

    •  liveset: The updated set of walkers.

    •  ns_params: The updated nested sampling parameters.

  ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  nested_sampling_loop!(liveset::AtomWalkers, n_steps::Int64, mc_routine::MixedMCRoutine, save_strategy::DataSavingStrategy)

  Perform a nested sampling loop for a given number of steps.

  Arguments
  ≡≡≡≡≡≡≡≡≡

    •  liveset::AtomWalkers: The initial set of walkers.

    •  n_steps::Int64: The number of steps to perform.

    •  mc_routine::MixedMCRoutine: The mixed Monte Carlo routine to use.

    •  save_strategy::DataSavingStrategy: The strategy for saving data.

  Returns
  ≡≡≡≡≡≡≡

    •  df::DataFrame: The data frame containing the iteration number and maximum energy for each step.

    •  liveset::AtomWalkers: The updated set of walkers.

    •  mc_routine.ns_params_main: The updated nested sampling parameters for the main routine.

Accessing the documentation

help?> NestedSamplingParameters
search: NestedSamplingParameters LatticeNestedSamplingParameters nested_sampling_step!

  mutable struct NestedSamplingParameters <: SamplingParameters

  The NestedSamplingParameters struct represents the parameters used in the nested sampling scheme.

  Fields
  ≡≡≡≡≡≡

    •  mc_steps::Int64: The number of total Monte Carlo moves to perform.

    •  initial_step_size::Float64: The initial step size, which is the fallback step size if MC routine fails to accept a move.

    •  step_size::Float64: The on-the-fly step size used in the sampling process.

    •  step_size_lo::Float64: The lower bound of the step size.

    •  step_size_up::Float64: The upper bound of the step size.

    •  fail_count::Int64: The number of failed MC moves in a row.

    •  allowed_fail_count::Int64: The maximum number of failed MC moves allowed before resetting the step size.

Finding phase transitions

# Calculate the ω factors
ωi = ωᵢ(energies.iter, 120)
# Shift the energies to be greater than or equal to zero
Ei = energies.emax .- minimum(energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:1000)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 3×6 for the 6-particle system
dof = 18
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the internal energies for each temperature
U = [internal_energy(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
cvs = cv(energies, β, dof, 120)
9991-element Vector{Float64}:
 0.0013175732622637696
 0.0013163656401052968
 0.0013164865839072522
 0.0013180887545059053
 0.0013209878572451813
 0.0013247504441241784
 0.0013288500299832545
 0.0013327866085318337
 0.0013361503619112659
 0.0013386436487862145
 0.0013400792149547592
 0.001340367194218547
 0.001339497919186227
 ⋮
 0.000865730605363293
 0.0008657001748975948
 0.0008656697600629232
 0.0008656393608486451
 0.0008656089772439815
 0.0008655786092383255
 0.0008655482568209475
 0.000865517919981241
 0.0008654875987085314
 0.0008654572929921685
 0.0008654270028215192
 0.0008653967281859931

Finding surface phase transitions

Finding surface phase transitions

Surface stuff

using AtomsBase
configs =  read_configs("slab.extxyz",pbc="TTF")
configs = FastSystem.(configs)
slabs = AtomWalker{2}.(configs, list_num_par=[80,6],frozen=[1,0])
480-element Vector{AtomWalker{2}}:
 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 ⋮
 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

 AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : 0.0 eV
    iter               : 0
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : 0.0 eV)

Surface stuff

slabs[1].configuration
FastSystem(H₈₆, periodic = TTF):
    bounding_box      : [ 11.2246        0        0;
                                0  9.72081        0;
                                0        0  29.1649]u"Å"

          .---------------------------.  
         /|                           |  
        / |                           |  
       /  |                           |  
      /   |                           |  
     /    |                           |  
    *     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     | H        H                |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |         H         H       |  
    |     |                H          |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     | H      H      H      H    |  
    |     |         H                 |  
    |  H  |  HH     HH     HH     H   |  
    |    H|     H      H      H       |  
    |   H |    H      H      H        |  
    H    HHH    HHH    HHH    HH      |  
    |     |                           |  
    |H    |HH     HH     HH     H     |  
    | H   |  H      H      H          |  
    |     | H      H      H      H    |  
    |  HH .--HHH----HHH----HHH----H---.  
    |    /                           /   
    |   HH     HH     HH     HH     /    
    H  /  HH     HH     HH     H   /     
    | /                           /      
    |H      H      H      H      /       
    *---------------------------*        

Surface stuff

surf_lj = LJParameters(epsilon=0.1, sigma=2.5, cutoff=4.0)
surf_ls = LJAtomWalkers(slabs, surf_lj)
mc = MCRandomWalkClone()
surf_ns_params = NestedSamplingParameters(500, 0.1, 0.01, 1e-5, 1.0, 0, 200)
save = SaveEveryN(n_traj=10, n_snap=50_000)
surf_energies, liveset, _ = nested_sampling_loop!(surf_ls, surf_ns_params, 50_000, mc, save)
(50000×2 DataFrame
   Row  iter   emax          Int64  Float64      
───────┼─────────────────────
     1 │     1    5.71085e12
     2 │     2    4.59391e9
     3 │     3    4.44588e9
     4 │     4    7.43776e8
     5 │     5    4.21399e8
     6 │     6    8.31856e7
     7 │     7    2.29957e7
     8 │     8    2.0287e7
     9 │     9    1.1382e7
    10 │    10    8.46159e6
    11 │    11    7.13711e6
   ⋮   │   ⋮         ⋮
 49991 │ 49991  -58.5274
 49992 │ 49992  -58.5274
 49993 │ 49993  -58.5274
 49994 │ 49994  -58.5274
 49995 │ 49995  -58.5274
 49996 │ 49996  -58.5274
 49997 │ 49997  -58.5274
 49998 │ 49998  -58.5274
 49999 │ 49999  -58.5274
 50000 │ 50000  -58.5274
           49979 rows omitted, LJAtomWalkers(AtomWalker{2}, LJParameters):
[1] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52736668197494 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[2] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.527366976199524 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[3] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52736700308327 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[4] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52736700465829 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[5] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.527367183158844 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
⋮
Omitted 470 walkers
⋮
[476] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.527535068052806 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[477] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52753654030334 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[478] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52753791580169 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[479] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.527562374940956 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
[480] AtomWalker{2}(
    configuration      : FastSystem(H₈₆, periodic = TTF, bounding_box = [[11.224620483093732, 0.0, 0.0], [0.0, 9.72080648619833, 0.0], [0.0, 0.0, 29.164864246657352]]u"Å")
    energy             : -58.52744314188147 eV
    iter               : 50000
    list_num_par       : [80, 6]
    frozen             : Bool[1, 0]
    energy_frozen_part : -54.68859657133728 eV)
LJParameters(0.1 eV, 2.5 Å, 4.0, -9.763240814208984e-5 eV)
, NestedSamplingParameters(500, 0.1, 0.009679022355186586, 1.0e-5, 1.0, 0, 200))

Surface results

surf_ls.walkers[1].configuration
FastSystem(H₈₆, periodic = TTF):
    bounding_box      : [ 11.2246        0        0;
                                0  9.72081        0;
                                0        0  29.1649]u"Å"

          .---------------------------.  
         /|                           |  
        / |                           |  
       /  |                           |  
      /   |                           |  
     /    |                           |  
    *     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |     |                           |  
    |    H|                           |  
    |     |H                    H     |  
    |     |                           |  
    |     | H      H      HH     HH   |  
    |     |                  H        |  
    |  H  |  HH     HH     HH     H   |  
    |    H|     H      H      H       |  
    |   H |    H      H      H        |  
    H    HHH    HHH    HHH    HH      |  
    |     |                           |  
    |H    |HH     HH     HH     H     |  
    | H   |  H      H      H          |  
    |     | H      H      H      H    |  
    |  HH .--HHH----HHH----HHH----H---.  
    |    /                           /   
    |   HH     HH     HH     HH     /    
    H  /  HH     HH     HH     H   /     
    | /                           /      
    |H      H      H      H      /       
    *---------------------------*        

Surface results

plot(surf_energies.iter[1000:end], surf_energies.emax[1000:end], ylabel="Energy", xlabel="Iteration", label="")

Surface results

# Calculate the ω factors
ωi = ωᵢ(surf_energies.iter, 480)
# Shift the energies to be greater than or equal to zero
Ei = surf_energies.emax .- minimum(surf_energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:1000)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 3×6 for the 6-particle system
dof = 18
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
surf_cvs = cv(surf_energies, β, dof, 480)
9991-element Vector{Float64}:
 0.0013842345024584256
 0.0014037564760710597
 0.001422094212288646
 0.0014393991765428245
 0.0014549553127525684
 0.0014679912279158577
 0.0014781013099614134
 0.0014853322720637684
 0.0014900807662374433
 0.0014929308416672571
 0.0014945083183555742
 0.001495382151241112
 0.0014960136907975394
 ⋮
 0.0030352937093422747
 0.003035254740728882
 0.0030352155242118813
 0.003035176059862044
 0.003035136347750496
 0.003035096387948459
 0.0030350561805273834
 0.003035015725558813
 0.0030349750231146062
 0.0030349340732666397
 0.003034892876087085
 0.003034851431648211

Surface results

Lattice systems

using Distributions, DataFrames

square_supercell_dimensions = (4, 4, 1)
square_basis=[(0.0, 0.0, 0.0)]

adsorption_energy = -0.04
nn_energy = -0.01
nnn_energy = -0.0025

# Initialize the lattice
occupied_sites = sample(1:16, 6, replace=false)

initial_lattice = LatticeSystem{SquareLattice}(;
           supercell_dimensions = square_supercell_dimensions,
           occupations=occupied_sites)

h = GenericLatticeHamiltonian(adsorption_energy, [nn_energy, nnn_energy], u"eV")
GenericLatticeHamiltonian{2,Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(eV,), 𝐋² 𝐌 𝐓⁻², nothing}}}:
    on_site_interaction:      -0.04 eV
    nth_neighbor_interactions: [-0.01, -0.0025] eV

Lattice systems

walkers = [deepcopy(initial_lattice) for i in 1:2000]

for walker in walkers
    walker.occupations = [false for i in 1:length(walker.occupations)]
    for i in sample(1:length(walker.occupations), 6, replace=false)
        walker.occupations[i] = true
    end
end

unique!(x -> x.occupations, walkers) # remove duplicates

# Now, we can create the liveset by combining the walkers and the potential
ls = LatticeGasWalkers(deepcopy(LatticeWalker.(walkers[1:1000])), h, perturb_energy=1e-10)
# We can now create a Monte Carlo object
mc = MCNewSample()
# Let's specify the nested sampling parameters: 
# Monte Carlo steps, energy perturbation, current fail count, allowed fail count
# Use `?LatticeNestedSamplingParameters` to see the documentation for more information
ns_params = LatticeNestedSamplingParameters(1,1e-10,0,500_000)
# We can also create a save object, using default value for most parameters except `n_snap` for saving snapshots
save = SaveEveryN(n_traj=10000, n_snap=500_000)

# And finally, we can run the simulation
lattice_energies, ls, _ = nested_sampling_loop!(ls, ns_params, 10_000, mc, save; accept_same_config=true)

Lattice systems

ωi = ωᵢ(lattice_energies.iter, 1000)
# Shift the energies to be greater than or equal to zero
Ei = lattice_energies.emax .- minimum(lattice_energies.emax)
# Specify the temperatures that we are interested in
Ts = collect(1:0.1:200)
# Define the Boltzmann constant
kb = 8.617333262e-5 # eV/K
# Calculate the inverse temperatures
β = 1 ./(kb.*Ts)
# Define the degrees of freedom, which is 2×4 for the 4-particle system
dof = 8
# Calculate the partition functions for each temperature
Zs = [partition_function(b, ωi, Ei) for b in β]
# Calculate the heat capacities as a function of temperature
lattice_cvs = cv(lattice_energies, β, dof, 1000)
1991-element Vector{Float64}:
 0.00034469333056180083
 0.00034469333142480505
 0.0003446933376293437
 0.00034469336959996606
 0.0003446934965557751
 0.0003446939063884268
 0.0003446950258828272
 0.0003446976939299924
 0.0003447033750094083
 0.00034471438632263944
 0.00034473410463458824
 0.00034476711998688886
 0.0003448193116499381
 ⋮
 0.00036374954549554323
 0.0003637334434403411
 0.0003637173611604454
 0.0003637012986247458
 0.00036368525580218975
 0.0003636692326617828
 0.0003636532291725878
 0.00036363724530372545
 0.00036362128102437375
 0.00036360533630376795
 0.00036358941111120074
 0.0003635735054160216

Lattice results

What’s next

  • Publish! There are some journals that are specifically designed for software papers

Journal of Open Source Software

Or, traditional journals like:

The Journal of Chemical Physics

Take-home message

  • Version control is a powerful tool
  • It’s a good practice for code development
  • It’s essential for collaborative projects
  • It has great importance in scientific reproducibility

Start using it, now!